Ειδικά Μαθήματα Βοτανικής

logo

Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε στον ιστότοπο του μαθήματος στο e-class

1 Πρώτα βήματα

1.1 Εγκατάσταση βιβλιοθήκης

Εγκατάσταση του πρώτου και κύριου πακέτου (και ενός βοηθητικού)
(Αυτό γίνεται μια φορά σε κάθε Η/Υ - εάν έχει εγκατασταθεί δεν χρειάζεται να δοθεί ξανά η εντολή αυτή)

## CRAN repository:
install.packages("easypackages", repos = "http://cran.us.r-project.org")

1.2 Φόρτωση βιβλιοθήκης

Φόρτωση του εν λόγω πακέτου

library(easypackages)

1.3 Εγκατάσταση βιβλιοθηκών

Με την εντολή packages('insert_package_name_here') γίνεται η εγκατάσταση των πακέτων που μας ενδιαφέρουν

Δοκιμάστε το

1.4 Φόρτωση βιβλιοθηκών

Εντολή για φόρτωση των πακέτων

libraries("rio", "car", "MASS", "psych", "leaps", "ggplot2", "klaR", "DiscriMiner", 
    "dplyr", "car", "MASS", "corrgram", "tidyverse", "plyr", "ggstatsplot")

1.5 Working Directory

Με την ακόλουθη εντολή, βλέπουμε σε ποιόν φάκελο δουλεύουμε:

getwd()

1.6 Αλλαγή του working directory

Ενώ με την εντολή αυτή, αλλάζουμε τον φάκελο στον οποίο δουλεύουμε και αποθηκεύουμε δεδομένα
[π.χ., setwd("E:/") εάν θέλουμε να δουλεύουμε στον σκληρό δίσκο Ε]

setwd("E:/Academic/Lessons/Labs/Labs")

1.7 Εισαγωγή αρχείου

Τώρα ας εισάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε

data(iris)

Όπως αντιλαμβάνεστε, σήμερα θα δουλέψουμε με δεδομένα τα οποία ήδη υπάρχουν εντός της R. Όμως για την εργασία που θα παραδώσετε εσείς, θα χρειαστεί να φορτώσετε δεδομένα από ένα αρχείο excel, το οποίο θα σας δοθεί στο τέλος της άσκησης. Για να το κάνετε αυτό, θα χρειαστεί να χρησιμοποιήσετε την βιβλιοθηκη readxl και την εντολη read_excel1

Ας δούμε τα δεδομένα μας:

names(iris)
## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width" 
## [5] "Species"
head(iris)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
5.1 3.5 1.4 0.2 setosa
4.9 3.0 1.4 0.2 setosa
4.7 3.2 1.3 0.2 setosa
4.6 3.1 1.5 0.2 setosa
5.0 3.6 1.4 0.2 setosa
5.4 3.9 1.7 0.4 setosa 
str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
iris$Species
##   [1] setosa     setosa     setosa     setosa     setosa     setosa    
##   [7] setosa     setosa     setosa     setosa     setosa     setosa    
##  [13] setosa     setosa     setosa     setosa     setosa     setosa    
##  [19] setosa     setosa     setosa     setosa     setosa     setosa    
##  [25] setosa     setosa     setosa     setosa     setosa     setosa    
##  [31] setosa     setosa     setosa     setosa     setosa     setosa    
##  [37] setosa     setosa     setosa     setosa     setosa     setosa    
##  [43] setosa     setosa     setosa     setosa     setosa     setosa    
##  [49] setosa     setosa     versicolor versicolor versicolor versicolor
##  [55] versicolor versicolor versicolor versicolor versicolor versicolor
##  [61] versicolor versicolor versicolor versicolor versicolor versicolor
##  [67] versicolor versicolor versicolor versicolor versicolor versicolor
##  [73] versicolor versicolor versicolor versicolor versicolor versicolor
##  [79] versicolor versicolor versicolor versicolor versicolor versicolor
##  [85] versicolor versicolor versicolor versicolor versicolor versicolor
##  [91] versicolor versicolor versicolor versicolor versicolor versicolor
##  [97] versicolor versicolor versicolor versicolor virginica  virginica 
## [103] virginica  virginica  virginica  virginica  virginica  virginica 
## [109] virginica  virginica  virginica  virginica  virginica  virginica 
## [115] virginica  virginica  virginica  virginica  virginica  virginica 
## [121] virginica  virginica  virginica  virginica  virginica  virginica 
## [127] virginica  virginica  virginica  virginica  virginica  virginica 
## [133] virginica  virginica  virginica  virginica  virginica  virginica 
## [139] virginica  virginica  virginica  virginica  virginica  virginica 
## [145] virginica  virginica  virginica  virginica  virginica  virginica 
## Levels: setosa versicolor virginica

Και ας μετονομάσουμε τα ονόματα των ειδών2:

iris$Species <- revalue(iris$Species, c(setosa = "Iris setosa", virginica = "Iris virginica", 
    versicolor = "Iris versicolor"))

1.8 Λιγότερα δεκαδικά

Με την ακόλουθη εντολή, θα εμφανίζονται πλέον λιγότερα δεκαδικά:

options(digits = 2)

1.9 Στατιστικά στοιχεία για κάθε taxon

Θέλουμε να μάθουμε την ελάχιστη, την μέγιστη και την μέση τιμή, καθώς και την τυπική απόκλιση3:

iris %>% split(.$Species) %>% map(~summary(.))
## $`Iris setosa`
##   Sepal.Length  Sepal.Width   Petal.Length   Petal.Width  
##  Min.   :4.3   Min.   :2.3   Min.   :1.00   Min.   :0.10  
##  1st Qu.:4.8   1st Qu.:3.2   1st Qu.:1.40   1st Qu.:0.20  
##  Median :5.0   Median :3.4   Median :1.50   Median :0.20  
##  Mean   :5.0   Mean   :3.4   Mean   :1.46   Mean   :0.25  
##  3rd Qu.:5.2   3rd Qu.:3.7   3rd Qu.:1.57   3rd Qu.:0.30  
##  Max.   :5.8   Max.   :4.4   Max.   :1.90   Max.   :0.60  
##             Species  
##  Iris setosa    :50  
##  Iris versicolor: 0  
##  Iris virginica : 0  
##                      
##                      
##                      
## 
## $`Iris versicolor`
##   Sepal.Length  Sepal.Width   Petal.Length  Petal.Width  
##  Min.   :4.9   Min.   :2.0   Min.   :3.0   Min.   :1.00  
##  1st Qu.:5.6   1st Qu.:2.5   1st Qu.:4.0   1st Qu.:1.20  
##  Median :5.9   Median :2.8   Median :4.3   Median :1.30  
##  Mean   :5.9   Mean   :2.8   Mean   :4.3   Mean   :1.33  
##  3rd Qu.:6.3   3rd Qu.:3.0   3rd Qu.:4.6   3rd Qu.:1.50  
##  Max.   :7.0   Max.   :3.4   Max.   :5.1   Max.   :1.80  
##             Species  
##  Iris setosa    : 0  
##  Iris versicolor:50  
##  Iris virginica : 0  
##                      
##                      
##                      
## 
## $`Iris virginica`
##   Sepal.Length  Sepal.Width   Petal.Length  Petal.Width  
##  Min.   :4.9   Min.   :2.2   Min.   :4.5   Min.   :1.40  
##  1st Qu.:6.2   1st Qu.:2.8   1st Qu.:5.1   1st Qu.:1.80  
##  Median :6.5   Median :3.0   Median :5.5   Median :2.00  
##  Mean   :6.6   Mean   :3.0   Mean   :5.6   Mean   :2.03  
##  3rd Qu.:6.9   3rd Qu.:3.2   3rd Qu.:5.9   3rd Qu.:2.30  
##  Max.   :7.9   Max.   :3.8   Max.   :6.9   Max.   :2.50  
##             Species  
##  Iris setosa    : 0  
##  Iris versicolor: 0  
##  Iris virginica :50  
##                      
##                      
## 
groupStds(iris[, 1:4], iris[, 5])
##              Iris setosa Iris versicolor Iris virginica
## Sepal.Length        0.35            0.52           0.64
## Sepal.Width         0.38            0.31           0.32
## Petal.Length        0.17            0.47           0.55
## Petal.Width         0.11            0.20           0.27

Μπορείτε επίσης να χρησιμοποιήσετε το σύνολο εντολών από την πρώτη εργαστηριακή άσκηση, προκειμένου να δημιουργήσετε τα στατιστικά αυτά δεδομένα.
Δοκιμάστε το για το IQR (Interquantile Range).

1.10 Γραφική απεικόνιση των αποτελεσμάτων

Ας αναπαριστήσουμε γραφικά τα δεδομένα μας για κάθε taxon και μια μεταβλητή4:

## Violin plot
ggplot(iris, aes(Species, Petal.Length)) + geom_violin(adjust = 4) + geom_boxplot(width = 0.1, 
    fill = "black", outlier.colour = NA) + stat_summary(fun.y = median, geom = "point", 
    fill = "white", shape = 21, size = 2.5) + xlab("Taxon") + ylab("Petal Length") + 
    theme(axis.title = element_text(face = "bold"), axis.text.x = element_text(face = "italic", 
        colour = "black"))

getwd()
## [1] "G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/EMB/2020"
ggsave("Petal length.png", width = 20, height = 20, units = "cm", dpi = 600)

Ένας άλλος τρόπος να οπτικοποιήσουμε τα δεδομένα μας είναι ο ακόλουθος5:

##--------------------------------------------------------------------------
## Another way to visualise the results
##--------------------------------------------------------------------------
ggstatsplot::ggbetweenstats(data = iris, x = Species, y = Petal.Length, mean.plotting = TRUE, 
    type = "robust", xlab = "Taxon", ylab = "Petal Length", title = "Dataset: Iris data set", 
    messages = FALSE)

Κάντε το ίδιο με τις υπόλοιπες μεταβλητές και αλλάξτε το όνομα του y άξονα.

ΠΡΟΣΟΧΗ
Για την εργασία σας, θα χρειαστεί να ελέγξετε εάν οι μεταβλητές σας έχουν κανονική κατανομή και εάν εμφανίζουν μεγάλο συντελεστή συσχέτισης. Ανατρέξτε σε προηγούμενες εργαστηριακές ασκήσεις για να δείτε πώς μπορείτε να το κάνετε αυτό.

2 Κυρίως ανάλυση

Είμαστε πλέον σε θέση να προχωρήσουμε στην κυρίως ανάλυση των μορφομετρικών μας δεδομένων. Είναι το πιο δύσκολο και κρίσιμο βήμα και θα χρειαστεί να είστε εξαιρετικά προσεκτικοί όταν θα τρέξετε την ακόλουθη εντολή για την εργασία σας:

attach(iris)  ## For your assignment, change that to mydata

iris <- iris[apply(iris,1,function(row) all(!is.na(row))),] ## Do the same here


set.seed(2934) ## Change that if you want to whichever number you like

cv.grp <- sample(1:10, nrow(iris), 
                 replace = TRUE) ## Remember to change the name

cv.lda <- c() ##A function that splits our data to 2 independent subsets, one to
              ## train our model and one to test how well our model behaves

for(i in 1:10) {
  train.ind <- which(cv.grp!=i)
  test.ind <- which(cv.grp==i)
  model.train <- lda(Species ~ Sepal.Length +## Change the variable names
                                             ## according to the ones that 
                                             ## appear in your data

  Sepal.Width + Petal.Width + Petal.Length,  ## Do the same here
  data = iris, subset = train.ind)           ## Change to mydata
  cv.lda[test.ind] <- predict(model.train,
  newdata = iris[test.ind,])$class             ## Do the same here
 }

mean(cv.lda==as.numeric(iris$Species))       ## And here
## [1] 0.98
Lec.lda <- lda(Species ~ Sepal.Length +      ## This is the test point
Sepal.Width + Petal.Width + Petal.Length,    ## Change the variable names

data = iris) ## Change to mydata


plot(Lec.lda,col = c("green","blue","red", "black", "orange") 
     
[iris$Species], pch = 19) ## Change to mydata

Με την ακόλουθη εντολή βλέπουμε τα αποτελέσματα της ανάλυσης μας:

Lec.lda


## Let's say what we did
##--------------------------------------------------------------------------------##
# Call: lda(Species ~ Sepal.Length + Sepal.Width + Petal.Width +
# Petal.Length, data = iris)
##--------------------------------------------------------------------------------##
## How much of the data lie in each species

## Prior probabilities of groups:

## setosa versicolor virginica

## 0.33 0.33 0.33
##--------------------------------------------------------------------------------##

## Group means:

## Sepal.Length Sepal.Width Petal.Width Petal.Length

## setosa 5.0 3.4 0.25 1.5

## versicolor 5.9 2.8 1.33 4.3

## virginica 6.6 3.0 2.03 5.6
##--------------------------------------------------------------------------------##

## Coefficients of linear discriminants:

## LD1 LD2

## Sepal.Length 0.83 0.024

## Sepal.Width 1.53 2.165

## Petal.Width -2.81 2.839

## Petal.Length -2.20 -0.932

##--------------------------------------------------------------------------------##

## Proportion of trace:

## The most important part

## It shows how many discriminant axes are there and which is the most
## important in order to distinguish the taxa

## The axes are like this: a*x1 + b*x2 We get a & b from the Coefficients of
## linear discriminants


## (here: 0.83*Sepal.Length + 1.53*Sepal.Width - 2.81*Petal.Width -
## 2.20*Petal.Length for the first axis)


## We can see them with: Lec.lda$scaling
##--------------------------------------------------------------------------------##

## LD1 LD2

## 0.9912 0.0088

## LD1 is the most important as it explains 99,12%
##--------------------------------------------------------------------------------##

2.1 Σημαντικότητα του κάθε διαγνωστικού άξονα σε σχέση με τους υπόλοιπους

Θέλουμε η αναλογία των τιμών αυτών, να είναι όσο το δυνατόν μεγαλύτερη, καθώς υποδηλώνουν τη σχετική σημαντικότητα κάθε άξονα. Όσο μεγαλύτερη η τιμή, τόσο πιο διαγνωστικός ο κάθε άξονας:

SVD <- Lec.lda$svd
SVD
## [1] 48.6  4.6

2.2 Συσχέτιση των μορφολογικών γνωρισμάτων με τον κάθε διαγνωστικό άξονα

Για την εργασία σας, εάν έχετε περισσότερους από 2 διαγνωστικούς άξονες, θα χρειαστεί να τους προσθέσετε:

LD1 <- predict(Lec.lda)$x[, 1]
LD2 <- predict(Lec.lda)$x[, 2]

## We will see which variable is the most important for each axis We are
## interested in the absolute value

cor(iris[, 1:4], LD1)  ## Change to mydata
##               [,1]
## Sepal.Length -0.79
## Sepal.Width   0.53
## Petal.Length -0.98
## Petal.Width  -0.97
## Change the number of columns

cor(iris[, 1:4], LD2)  ## Do the same
##               [,1]
## Sepal.Length 0.218
## Sepal.Width  0.758
## Petal.Length 0.046
## Petal.Width  0.223

2.3 Πίνακας Ορθής Ταξινόμησης

Σε αυτό το σημείο, θέλουμε να δούμε πόσο καλά διακρίνουν οι μορφολογικοί χαρακτήρες οι οποίοι αναδείχθηκαν ως σημαντικοί από την ανάλυση μας, τα διάφορα taxa μεταξύ τους:

x <- Lec.lda
y <- predict(x, iris)  ## Change to mydata

detach(iris)  ## Change to mydata - DON'T FORGET TO DO THAT

round(100 * errormatrix(iris$Species, y$class, relative = T), 1)  ## Change to mydata
##                  predicted
## true              Iris setosa Iris versicolor Iris virginica -SUM-
##   Iris setosa             100               0              0     0
##   Iris versicolor           0              96              4     4
##   Iris virginica            0               2             98     2
##   -SUM-                     0              33             67     2

2.4 Σημαντικότητα του κάθε μορφολογικού χαρακτήρα

Τώρα μπορούμε να διαπιστώσουμε πόσο σημαντικός είναι ο εκάστοτε μορφολογικός χαρακτήρας:

## Change to mydata

## The smaller the Wilk's lambda, the better

greedy.wilks(Species ~ Sepal.Length + Sepal.Width + Petal.Width + Petal.Length, 
    data = iris, niveau = 0.1)
## Formula containing included variables: 
## 
## Species ~ Petal.Length + Sepal.Width + Petal.Width + Sepal.Length
## <environment: 0x000000002143c1a0>
## 
## 
## Values calculated in each step of the selection procedure: 
## 
##           vars Wilks.lambda F.statistics.overall p.value.overall
## 1 Petal.Length        0.059                 1180         2.9e-91
## 2  Sepal.Width        0.037                  307        2.8e-103
## 3  Petal.Width        0.025                  258        5.0e-113
## 4 Sepal.Length        0.023                  199        1.4e-112
##   F.statistics.diff p.value.diff
## 1            1180.2      2.9e-91
## 2              43.0      2.0e-15
## 3              34.6      5.1e-13
## 4               4.7      1.0e-02

2.5 Γραφική απεικόνιση των αποτελεσμάτων

s <- predict(Lec.lda)$x
s <- as.data.frame(s)
Taxon <- factor(iris$Species)  ## ## Change to mydata 

s <- cbind(Taxon, s)
head(s, 1)
Taxon LD1 LD2
Iris setosa 8.1 0.3 
names(s)
## [1] "Taxon" "LD1"   "LD2"
getwd()
## [1] "G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/EMB/2020"
ggplot(s, aes(x = LD1)) + geom_histogram(fill = "white", colour = "black") + 
    facet_grid(Taxon ~ .)

ggsave("Iris LD1.png", width = 20, height = 20, units = "cm", dpi = 600)

ggplot(s, aes(x = LD2)) + geom_histogram(fill = "white", colour = "black") + 
    facet_grid(Taxon ~ .)

ggsave("Iris LD2.png", width = 20, height = 20, units = "cm", dpi = 600)


a <- s$LD1
b <- s$LD2
TAXON <- iris$Species  ## ## Change to mydata

g <- data.frame(TAXON, a, b)  ## Create a dataframe that contains all the above

## Another dataframe containing the centroids for each species and each
## discriminant axis
centroids <- aggregate(cbind(a, b) ~ TAXON, g, mean)

ggplot(g, aes(a, b, color = TAXON, shape = TAXON)) + geom_point(size = 3) + 
    stat_ellipse(geom = "polygon", alpha = 1/2, aes(fill = TAXON)) + xlab("LD1") + 
    ylab("LD2") + theme(axis.title = element_text(face = "bold")) + stat_summary(data = centroids, 
    fun.y = "mean", geom = "point", shape = 22, size = 4, colour = "black", 
    aes(fill = TAXON)) + 
## Shape of each species

## If you have more than 3 species you must add some values (range: 10-20)
scale_shape_manual(values = c(15, 16, 17))

ggsave("Iris ellipses.png", width = 20, height = 20, units = "cm", dpi = 600)

3 Εργασία για το σπίτι

Έχετε διορία επτά (7) ημερών να στείλετε την εργασία σας σε μορφή PDF6 σε αυτό το e-mail και να απαντήσετε στα ακόλουθα ερωτήματα, αφού έχετε περιγράψει το σετ δεδομενων το οποίο έχετε χρησιμοποιήσει:
1. Πόσα είδη έχει;
2. Πόσα δειγματα έχει το κάθε είδος;
3. Πόσους μορφολογικούς χαρακτήρες περιέχει το σετ δεδομένων σας;
4. Με βάση τα αποτελέσματα σας, να φτιάξετε έναν πίνακα ο οποίος να περιέχει τη μέση, την ελάχιστη και την μέγιστη τιμή για κάθε μορφολογικό χαρακτήρα κάθε είδους, μαζί με την τυπική απόκλιση.
5. Ποιές μεταβλητές δεν εμφανίζουν κανονική κατανομή;
6. Ποιές μεταβλητές εμφανίζουν μεγάλο (\(>0.8\)) συντελεστή συσχέτισης;
7. Εάν υπάρχουν κάποια ζεύγη ανεξάρτητων μεταβλητών με σ.σ. \(>0.8\), θα τις χρησιμοποιήσετε όλες;
8. Πόσοι διαγνωστικοί άξονες εμφανίστηκαν στην ανάλυση σας; Ποιό είναι το ποσοστό της ποικιλότητας που εξηγεί ο καθένας εξ αυτών; Ποιός είναι ο σημαντικότερος και πώς μπορείτε να είστε βέβαιοι για την απάντηση σας;
9. Ποιός μορφολογικός χαρακτήρας σχετίζεται περισσότερο με κάθε εξ αυτών των αξόνων;
10. Δημιουργήστε τον πίνακα ορθής ταξινόμησης. Ποιό είναι το ποσοστό ορθής ταξινόμησης για κάθε taxon το οποίο έχετε στην ανάλυση σας; Υπάρχει ποσοστό αλληλεπικάλυψης μεταξύ κάποιων taxa; Εάν ναι, που θεωρείτε ότι οφείλεται;
11. Ποιός μορφολογικός χαρακτήρας είναι ο σημαντικότερος για την διάκριση των taxa που έχετε στην ανάλυση σας;
12. Δημιουργήστε ένα γράφημα το οποίο να απεικονίζει γραφικά τα αποτελέσματα της ανάλυσης σας. Τι παρατηρείτε; Συμφωνούν τα αποτελέσματα του πίνακα ορθής ταξινόμησης με το εν λόγω γράφημα;

  1. Hint: Ανατρέξτε σε προηγούμενες εργαστηριακές ασκήσεις για να δείτε πώς γίνεται ή διαβάστε το manual της εν λόγω βιβλιοθήκης: readxl↩

  2. Δεν θα χρειαστεί να το κάνετε αυτό εσείς για την εργασία σας↩

  3. Για την εργασία σας, για να βρείτε την τυπική απόκλιση, θα χρειαστεί να τρέξετε την ακόλουθη εντολή πρώτα: mydata.sd <- as.data.frame(mydata) - εννοείται στην περίπτωση που έχετε ονομάσει το αρχείο σας mydata↩

  4. Για να κάνετε αυτό το βήμα για την εργασία σας, θα πρέπει να αλλάξετε το όνομα του αρχείου, καθώς και τα ονόματα των μεταβλητών - Εάν έχετε απορίες, τώρα είναι η στιγμή να τις εκφράσετε↩

  5. Για την εργασία σας μπορείτε να παρουσιάσετε όποιον από τους 2 τρόπους επιθυμείτε↩

  6. Διαφορετικά δεν θα γίνει δεκτή και δεν θα βαθμολογηθείτε↩

Κώστας Κουγιουμουτζής

Spring semester 2019-2020

 

Delivered to you by Dr. Kostas Kougioumoutzis

kkougiou@aua.gr